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t-H [ Abstract. Substantial amount of information both on the source and on character- 

f^ ' istics of intercepting magnetic fields is encoded in the distribution in arrival times, 

directions, and energies of charged ultra-high energy cosmic rays from discrete sources. 

We present a numerical approach that allows to extract such information from data 

from next generation experiments. 
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INTRODUCTION 



^ ■ The origin of ultra-high energy cosmic rays is still a major unresolved mystery in 

->^ . astrophysics. It is hard to imagine a mechanism producing particles of energy up 

'"5 I to several 100 EeV (= 10^° eV). In addition, sources must be closer than ^ 50Mpc 

I ' because of the limited range of nucleons due to photo-pion production at these en- 

;_i ■ ergies. No obvious astrophysical sources have been found within this distance [1,2]. 

c/3 . Although deflection of charged primaries can be strong in the direction of strong 

' magnetic fields along large mass agglomerations such as the supergalactic plane [3] , 

a deflection of several degrees at the most is expected along most other lines of sight 

^ ■ for nucleons above 10^° eV, due to the Faraday rotation limit on the large-scale field. 

For next generation experiments with their anticipated much improved exposure, 

this opens up the possibility to do "particle astronomy" and pinpoint sources along 

the arrival directions. 

It has been noted in that respect that a sub class of events above 4 x 10^^ eV 
seems to cluster in arrival directions [4]. If these clusters originated in discrete 
sources, some interesting qualitative consequences result already, such as a limit on 
the intercepted magnetic fields that is comparable to the Faraday rotation limit [5] . 
Next generation experiments should in this case see clusters of several tens or even 
hundreds of events at these energies in case of the Pierre Auger Project [6] and the 
Orbital Wide-angle Light Collector (OWL) [8], respectively. This just follows from 
scaling to the relevant expected exposures. A data pool of arrival directions, times, 
and energies of that size contains a substantial amount of information on both the 
source of a given cluster of events and magnetic fields intercepting the line of sight. 



This motivated us to conduct a detailed feasibility study for the potential of 
future experiments to reconstruct certain parameters that characterize the source 
mechanism and the large-scale magnetic field both of which are poorly known at 
present. We first briefly describe our method and then summarize our results and 
give some examples. 

DESCRIPTION OF APPROACH 

Here we describe the essential ingredients of our numerical approach; more details 
can be found in Refs. [9,10]. 

The propagation of nucleons through extragalactic space is simulated using the 
Monte Carlo technique: First, a magnetic field realization is set up on a grid via Fast 
Fourier Transformation by sampling a power spectrum of the form {B'^{k)) oc A;"^ 
for wavenumbers k < 2ti /l^. and otherwise, where Ic characterizes the coherence 
scale and ub the magnetic field power spectrum. As results are quite insensitive to 
ns, we assume n^ = if not stated otherwise. More sophisticated models for the 
magnetic field including the role of the large scale structure of galaxies, for example, 
along the lines discussed in Ref. [3], may be implemented in the future. Many 
nucleon trajectories are then calculated between a given source and observer by 
sampling direction of emission, injection energy and the stochastic pion production 
loss that becomes important above the Greisen-Zatsepin-Kuzmin (GZK) cut-off [11] 
at a few 10^^ eV . Pair production by protons has been incorporated as a continuous 
energy loss. One of the main problems that has to be solved when images of discrete 
sources are discussed, has not been considered in other work on propagation and 
deflection [12] and consists of the fact that different trajectories not only originate 
at the same source, but also have to reach the same observer. 

From the injection energies, direction, time, and energy of arrival recorded for 
the trajectories we then calculate histograms for the distribution in these quantities 
by convolving with the injection spectrum (typically a power law with index 7 for 
the differential spectrum) in energy and with a timescale To, that characterizes the 
emission timescale. Histograms are also smeared out in energy to account for finite 
energy resolution (typically ^E / E ~ 0.14, a value expected for future detectors) 
and are proportional to the source fiuence A^o- We also use the parameter 



rE - 2.0 {^\ i^—Y i^^] f^ yr. (1) 

^ l^SOMpcyl llOOEeVy Uo-ncy l^lMpc^^ ^' 

which is the average time delay for a proton of energy E over a distance D in a 
field of r.m.s. strength -Brms when energy loss is negligible, and te ^ D [13]. It is 
related to the average deflection angle 6e by 

/ n \ "^/^ / T \ ^''^ 



The subscript E is given in EeV in the following. 

Clusters of events are then obtained by sampling the histogram with Poisson 
statistics over a time window of width Tobs which constitutes the experimental 
lifetime, at a random position. Conversely, for a given cluster of events, a likeli- 
hood can be calculated for a given histogram that corresponds to certain values 
of the physical parameters described above. Averaging over different observational 
window positions and different realizations of the magnetic field for the same pa- 
rameters yields the likelihood function £(rioo, ^s, -D, 7, Nq, Ic, ub)- Marginalization 
over part of these parameters, using priors that account for certain constraints and 
other available information, can be used to reduce the parameter space. 



RESULTS AND EXAMPLES 

We first give a brief outline of the main features of the angle-time-energy images 
of clusters of ultra-high energy nucleons which have been described in detail in 
Ref. [14]. 

If both Ts < Tioo, and rioo is small compared to Tobs, arrival time and energy 
are correlated according to te oc E^'^; see Eq. (1). The angular image can not be 
resolved in this case. 

A source, such that rioo ^ Ts and rioo ^ Tobs, can be seen only in a limited range 
of energies, at a given time, as first pointed out in Ref. [13], and demonstrated in 
Fig. 1. Below the GZK cut-off, the width of this stripe, in the time-energy plane 
and within the observational window of length Tobs, is then governed by the ratio 
DOe/Ic for the energy at which events are observed: If this ratio is much smaller 
than 1, all nucleons have experienced the same magnetic field structure during their 
propagation and the width is very small in the absence of pion production; in the 
opposite case the width is expected to be Ate/te ~ 60%, even for negligible energy 
loss. Furthermore, the angular image is point-like or diffuse, with 6e describing the 
systematic off-set from the direction to the source, or the angular extent of the 
diffuse image that is centered on the source, respectively, in these two cases. If 
D9e/Ic ~ 1, several images of the source can result [15]. 

For a source emitting continuously at all energies of interest here, i.e. with 
^s ^ T30 and Ts ^ Tobs, events of any energy can be recorded at any time. 
Whereas the above remarks on the angular image now apply for all energies (see 
Fig. 2), the distribution of arrival time vs. energy is now uniform. 

Finally, for a source, such that rioo < ^s and r^o > Ts, together with Ts ^ Tobs, 
there exists an energy Eq, such that te^. = Ts. In this case, protons with E < Eq 
are not detected, as they could not have reached us within Tobs- However, protons 
with E > Eq are detected as for a continuously emitting source, i. e. with a uniform 
distribution of arrival times vs. energy (see, e.g.. Fig. 3). 

We now summarize results on the potential to reconstruct the parameters rioo, 
Ts, -D, 7, A^o, ^c, and n^ in these scenarios. Details have been presented in Ref. [10]. 



The likelihood presents different degeneracies between different parameters, 
which complicates the analysis. As an example, the likelihood is degenerate in 
the ratios Nq/T^, or Nq/Atioq, with A^o the total ffuence, and Arioo the spread in 
arrival time: these ratios represent rates of detection. Another example is given by 
the degeneracy between the distance D and the injection energy spectrum index 7. 
Yet another is the ratio D9e/Ic oc {Dte)^^'^/Ic, that controls the size of the scatter 
around the mean of the te — E correlation. Therefore, in most general cases, values 
for the different parameters cannot be pinned down, and generally, only domains of 
validity are found. We remark, however, that the generic scenarios discussed above 
are, in general, easy to distinguish from the likelihood function (see, e.g.. Fig. 2). 

We find that the distance to the source is obtained from the pion production sig- 
nature, above the GZK cut-off, when the emission timescale of the source dominates 
over the time delay. The lower the minimal energy above which the source appears 
as emitting continuously, the higher the accuracy on the distance D. The error on 
D is, in the best case, typically a factor 2, for one cluster of ~ 40 events. In this 
case, where the emission timescale dominates over the time delay at all observable 
energies, information on the magnetic field is only contained in the angular im- 
age. Qualitatively, the size of the angular image is proportional to -Brms(-D/c)^^^/-S, 
whereas the structure of the image, i.e. the number of separate images, is con- 
trolled by the ratio D9e/Ic oc D^^'^ B mis /Ell^'^. Finally, the case where the time 
delay dominates over the emission timescale, with a time delay shorter than the 
lifetime of the experiment, also allows to estimate the distance with a reasonable 
accuracy. 

The injection spectrum index 7 can be measured provided ultra-high energy 
cosmic rays are recorded over a bandpass in energy that is sufficiently broad. In 
general, it is comparably easy to rule out a hard injection spectrum if the actual 
7 ^ 2.0, but it is much harder to distinguish between 7 = 2.0 and 2.5. 

The strength of the magnetic field can only be obtained from the time-energy 
image in this latter case because the angular image will not be resolvable. When the 
time delay dominates over the emission timescale, and is, at the same time, larger 
than the lifetime Tobs of the experiment, only a lower limit corresponding to Tobs, 
can be placed on the time delay, hence on the strength of the magnetic field. When 
combined with the Faraday rotation upper limit, this would nonetheless allow to 
bracket the r.m.s. magnetic field strength within a few orders of magnitude. Here 
as well, significant information is contained in the angular image. 

The coherence length enters the ratio {DteY^'^/Ic that controls the scatter around 
the mean of the te — E correlation in the time-energy image. It can therefore 
be estimated from the width of this image, provided the emission timescale is 
dominated by te (otherwise the correlation would not be seen), and some prior 
information on D and te is available. If the source appears continuous and the 
time delay is large enough to resolve the angular image, l^ can be constrained or 
even estimated from the fact that DOe/Ic passes through 1 at the energy where 
the scatter AOe/Oe becomes comparable to 1 (it is much smaller at energies that 
are higher but still below the GZK cut-off; see Fig. 2). Our simulations showed no 



sensitivity to the magnetic field power spectrum characterized by hb- 

An emission timescale much larger than the experimental lifetime may be esti- 
mated if a lower cut-off in the spectrum is observable at an energy Eq, indicating 
that To, ~ teq- The latter may, in turn, be estimated from the angular image size 
via Eq. (2), where the distance can be estimated from the spectrum visible above 
the GZK cut-off, as discussed above. An example for this scenario is shown in 
Fig. 3. For angular resolutions A^, timescales in the range 

could be probed. The lower limit follows from the requirement that it should be 
possible to estimate Te from 6e, using Eq. (2), otherwise only an upper limit on 
Ts, corresponding to this same number, would apply. The upper bound in Eq. (3) 
comes from constraints on maximal time delays in cosmic magnetic fields, such 
as the Faraday rotation limit in case of the cosmological large-scale field (smaller 
number) and knowledge on stronger fields associated with the large-scale galaxy 
structure (larger number). Eq. (3) constitutes an interesting range of emission 
timescales for many conceivable scenarios of ultra-high energy cosmic rays. For 
example, the hot spots in certain powerful radio galaxies that have been suggested 
as ultra-high energy cosmic ray sources [16], have a size of only several kpc and 
could have an episodic activity on timescales of ~ 10^ yr. 

CONCLUSIONS 

A wealth of information on both the production mechanism of highest energy 
cosmic rays and on the structure of large-scale magnetic fields is encoded in angle- 
time-energy images of discrete sources. If the clustering suggested by AGASA is 
real, tens (for the Pierre Auger Project) to hundreds (for the OWL Project) of 
events above a few 10^^ eV can be expected from individual sources alone. With 
resolutions of 10-20% in energy and fractions of a degree in angle, next generation 
experiments should be able to exploit this information. 

Special thanks go to Martin Lemoine for an ongoing extensive collaboration on 
this subject. I also thank Angela Olinto and David Schramm for collaboration in 
earlier stages. This work was supported, in part, by the DoE, NSF, and NASA at 
the University of Chicago. 
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FIGURE 1. (a) An arrival timc-cncrgy histogram for 7 = 2.0, tiqo ~ lOOyr, Ts <C tioo, 
Zc — 1 Mpc, D = 50Mpc, corresponding to Brms — 4 x 10^^^ G. Contours are in steps of a factor 
lO^"* = 2.51; (b) Observable energy spectrum for several positions of the observational window in 
the histogram in (a); (c) Example for a cluster in the arrival time-energy plane resulting from one 
of the cuts shown in (b); (d) The likelihood function, marginalized over A^o and 7, for D = 50 Mpc, 
Ic = 0.25 Mpc, for the cluster shown in (c), in the Ts — tioo plane. The contours shown go from 
the maximum down to about 0.01 of the maximum in steps of a factor lO'^'^ = 1.58. The fall-off 
at Tioo ^ 50 yr and Ts ;$ Syr is a numerical artifact due to limited statistics. The true parameters 
are reasonably well reconstructed. 
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FIGURE 2. An angle-histogram for 7 = 2.0, noo = lO'^yr, Ts > noo, Ic ^ 1 Mpc, D = 50Mpc, 
corresponding to Brms — 4 x 10^^" G. An angular resolution of 0.05° was assumed. The point 
a; = y = is the source position and the contours decrease in steps of 0.1 in the logarithm to base 
10. (a) Image integrated over all energies E > 30EeV. Two partially blended, elongated images 
at x ~ 0.2°, y ~ 0.2° and at a; ~ 0.7°, y ~ 0.4° are clearly visible, the second one being more 
luminous by about a factor 4; (b) Same for E > 100 EeV. The two images are now closer to the 
source position; (c) Same for E > 200 EcV. The second image has almost disappeared; (d) Same 
for E > 300 EeV. As a consequence, D6e/Ic ~ 1 for i? ~ lOOEeV. If D can be estimated from 
the energy spectrum, an estimate for Ic results. 
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FIGURE 3. (a) An arrival time-energy histogram for 7 = 2.0, tioo = 50 yr, Tq — 200 yr, 
Zc — 1 Mpc, D = 50Mpc, corresponding to Brms — 3 x 10^^^ G. Contours are in steps of a factor 
lO*^-* = 2.51; (b) Example for a cluster in the arrival time-energy plane resulting from the cut 
indicated in (a) by the dashed line at t ~ 100 yr; (c) Same as Fig. 1 (d), but for the cluster shown 
in (b). Note that the likelihood clearly favors Ts ~ t^q. For tiqo large enough to be estimated 
from the angular image size, Ts S> Tobs can, therefore, be estimated as well. 



